========================================================================= 
Log Path: ./log/SpainReplicationLog.log 
Program Path: /Users/hannah/Dropbox/Research/Under Review/Local Suffrage/JOP/Replication/SpainReplication.R 
Working Directory: /Users/hannah/Dropbox/Research/Under Review/Local Suffrage/JOP/Replication 
User Name: hannah 
R Version: 4.3.0 (2023-04-21) 
Machine: ls-halarian-2019-iMac.local x86_64 
Operating System: Darwin 22.3.0 Darwin Kernel Version 22.3.0: Mon Jan 30 20:42:11 PST 2023; root:xnu-8792.81.3~2/RELEASE_X86_64 
Base Packages: parallel stats graphics grDevices utils datasets methods base 
Other Packages: tidylog_1.0.2 tjbal_0.4.0 kbal_0.1 kableExtra_1.3.4 modelsummary_1.4.1 fixest_0.11.1 readstata13_0.10.1 devtools_2.4.5 usethis_2.2.2 panelView_1.1.17 doParallel_1.0.17 iterators_1.0.14 foreach_1.5.2 lubridate_1.9.2 forcats_1.0.0 stringr_1.5.0 dplyr_1.1.4 purrr_1.0.1 readr_2.1.4 tidyr_1.3.0 tibble_3.2.1 ggplot2_3.4.4 tidyverse_2.0.0 logr_1.3.6 
Log Start Time: 2024-03-16 13:57:24.190795 
========================================================================= 

> #' ---
> #' title: "Local Suffrage Increases Citizenship Acquisition: Evidence from the European Union 
> #' Replication Log: Spain - Trajectory Balancing & Event Study"
> #' author: "Hannah M. Alarian"
> #' date : "March 16, 2024"
> #' ---
> 
> #Use file Spain.dta
> #Variable information can be found in Codebook tab Spain
> 
> #####Libraries-------------------------------------------------------------
> library(tidyverse)
> library(ggplot2)
> library(foreach)
> library(doParallel)
> require(parallel)
> library(panelView)
> library(devtools)
> library(readstata13)
> library(dplyr)
> library(fixest)
> library(modelsummary)
> library(kableExtra)
> 
> library(kbal) #install via github #install_github("chadhazlett/KBAL")
> library(tjbal) #install via github #install_github("xuyiqing/tjbal")
> 
> #Data Preparation--------------------------------------------------------------
> #Set working directory to location of data. Call in data from working directory
> #setwd()
> spain <- read.dta13("~/Dropbox/Research/Under Review/Local Suffrage/JOP/Data/Spain/Replication/Spain.dta")
> head(spain)
> 
> #Removing Missing Cases and Unbalanced
> spain <- spain[spain$scou != "BIH" & spain$scou != "BLR" & spain$scou != "ISR" &
>                  spain$scou != "SYR" & spain$scou != "IRQ" & spain$scou != "CUB" &
>                  spain$scou != "TUR" & spain$scou!= "GEO", ]
> 
> #Re-Coding Treatment variable 
> spain$treat <- 0
> spain$treat[spain$municipal == 1 & spain$year >= 2010] <- 1
> spain$treat[spain$scou =="CPV" & spain$year >= 2010] <- 1 #to run treatment at same time
> 
> #Classifying for Event Study
> spain$treated <- 0
> spain$treated[spain$scou == "ECU" | spain$scou == "BOL" | spain$scou == "PER" 
>             | spain$scou == "COL" | spain$scou ==  "CHL" | spain$scou =="PRY" 
>             | spain$scou == "CPV"] <- 1
> spainT <- spain[spain$treated == 1,]
> 
> #Relabeling variables for graphs
> spain <- spain %>%
>   rename("Acquisitions" = "citacq",
>     Residents = lpermres,
>     GDPpc = sgdp,
>     Unemployment = sunemploy,
>     Democracy = spolity2,
>     Colony = colony,
>     Language = comlang_ethno,
>     Perm.Res= permres
>   )
> 
> #Subsetting for balance prior to 2012 & up to 2015 & including 2007
> spain12 <- spain[spain$year >= 2008 & spain$year <= 2012, ]
> spain12 <- spain12[spain12$scou %in% names(which(table(spain12$scou) > 4)), ]
> 
> spain15 <- spain[spain$year >= 2008, ]
> spain15 <- spain15[spain15$scou %in% names(which(table(spain15$scou) > 7)), ]
> 
> spain7 <- spain[spain$year >= 2007, ]
> spain7 <- spain7[spain7$scou %in% names(which(table(spain7$scou) > 8)), ]
> 
> spainT <- spainT[spainT$year >= 2008 & spainT$year <= 2012, ]
> 
> #Looking at data structure
> panelview(Acquisitions ~ treat, data = spain12, show.id = c(1:49),
>           index = c("scou", "year"), xlab = "Year", ylab = "Sending Country",
>           axis.lab.gap = c(0,1), pre.post = TRUE)
> with(spainT, table(year, treat))
> 
> #Replicating Analysis--------------------------------------------------------------
> ##Models##
> #Average Treatment Effect on the Treated by Year (Mean Balanced)
> out.mbal <- tjbal(Y= "Acquisitions", D = "treat", 
>                   X = c("Colony", "Language","Unemployment"),
>                   data = spain12,
>                   X.avg.time = list(c(2008), c(2008), c(2008:2009)),
>                   estimator = "mean",
>                   index = c("scou", "year"), parallel = F,
>                   demean = T, vce = "boot", nsims = 1000, seed = 1013, conf.lvl = .9)
> #Trajectory Balancing - Maximizing Covariate Inclusion
> out.mbal1r <- tjbal(Y= "Acquisitions", D = "treat", 
>                     X = c("Colony", "Language","Unemployment", "Democracy", "GDPpc"),
>                     data = spain12,
>                     X.avg.time = list(c(2008), c(2008), c(2008:2009), c(2008:2009), c(2008:2009)), 
>                     estimator = "mean",
>                     index = c("scou", "year"), parallel = F,
>                     demean = T, vce = "boot", nsims = 1000, seed = 1013, conf.lvl = .9)
> #Trajectory Balancing with All Available Time Points (2007 - 2015)
> out.mbal7 <- tjbal(Y= "Acquisitions", D = "treat", 
>                    X = c("Colony", "Language","Unemployment", "Democracy", "GDPpc"),
>                    data = spain7,
>                    X.avg.time = list(c(2008), c(2008), c(2007:2009), c(2007:2009), c(2007:2009)), 
>                    estimator = "mean",
>                    index = c("scou", "year"), parallel = F,
>                    demean = T, vce = "boot", nsims = 1000, seed = 1013, conf.lvl = .9)
> #Permanent Residents Trajectory Balancing (2008 - 2015)
> out.mbal15r2 <- tjbal(Y= "Acquisitions", D = "treat", 
>                       X = c("Perm.Res"),
>                       data = spain15,
>                       X.avg.time = list(c(2008:2009)), 
>                       estimator = "mean",
>                       index = c("scou", "year"), parallel = F,
>                       demean = T, vce = "boot", nsims = 1000, seed = 1013, conf.lvl = .9)
> #Event Study
> set.seed(1013) #setting seed for replication
> M <- feols(citacq ~ i(year, ref = 2009), data = spainT,
>            cluster = 'scode')
> names(M$coefficients) <- c("2009", "2008", "2010", "2011", "2012")
> 
> #Figure 3--------------------------------------------------------------
> #Panel a
> plot(out.mbal, type = "ct", ylim = c(0, 15000), 
>      main = "Treated and Counterfactual Averages\n",
>      ylab = "Citizenship Acquisitions", xlab = "Year", count = F)
> #Panel b
> plot(out.mbal, ylab = "Citizenship Acquisitions", xlab = "Year", 
>      main = "Average treatment Effect on the Treated \n90% CI", count = F)
> 
> #Appendix B.1.2--------------------------------------------------------------
> print(out.mbal)
> 
> #Appendix C.1.2--------------------------------------------------------------
> #Panel a
> plot(out.mbal1r, type = "ct", ylim = c(0, 15000), main = "Treated and Counterfactual Averages\n",
>      ylab = "Citizenship Acquisitions", xlab = "Year", count = F)
> #Panel b
> plot(out.mbal1r, ylab = "Citizenship Acquisitions", xlab = "Year",
>      main = "Average treatment Effect on the Treated \n90% CI", count = F)
> #Table
> print(out.mbal1r)
> 
> #Appendix C.1.3--------------------------------------------------------------
> #Panel a
> plot(out.mbal7, type = "ct", ylim = c(0, 20000), main = "Treated and Counterfactual Averages\n",
>      ylab = "Citizenship Acquisitions", xlab = "Year", count = F)
> #Panel b
> plot(out.mbal7, ylab = "Citizenship Acquisitions", xlab = "Year",
>      main = "Average treatment Effect on the Treated \n90% CI", count = F)
> #Table
> print(out.mbal7)
> 
> #Appendix C.1.4--------------------------------------------------------------
> #Panel a
> plot(out.mbal15r2, type = "ct", ylim = c(0, 20000), main = "Treated and Counterfactual Averages\n",
>      ylab = "Citizenship Acquisitions", xlab = "Year", count = F)
> #Panel b
> plot(out.mbal15r2, ylab = "Citizenship Acquisitions", xlab = "Year",
>      main = "Average treatment Effect on the Treated \n90% CI", count = F)
> #Table
> print(out.mbal15r2)
> 
> #Appendix C.1.5--------------------------------------------------------------
> #Standardized Difference in Weighted and Unweighted Coefficients
> #Top Right panel
> plot(out.mbal, type = "balance", main = "Covariate Balance\nMean\n(2008-2012)")
> #Top Left panel
> plot(out.mbal1r, type = "balance", main = "Covariate Balance\nMean\n(2008-2012)")
> #Bottom Right
> plot(out.mbal7, type = "balance", main = "Covariate Balance\nMean\n(2007-2015)")
> #Bottom Left
> plot(out.mbal15r2, type = "balance", main = "Covariate Balance\nMean\n(2008-2015)")
> 
> 
> #Standardized Difference in Weighted and Unweighted Standard Deviations
> #Top Right panel
> sdm <- plot(out.mbal, type = "balance", stat = "sd",
>             main = "Covariate Balance\nMean\n(2008-2012)")
> sdm + scale_color_manual(values = c("lightcoral", "purple"))
> #Top Left panel
> sdmMax <- plot(out.mbal1r, type = "balance", stat = "sd",
>                main = "Covariate Balance\nMean\n(2008-2012)")
> sdmMax + scale_color_manual(values = c("lightcoral", "purple"))
> #Bottom Right
> sdk <- plot(out.mbal7, type = "balance", stat = "sd",
>             main = "Covariate Balance\nMean\n(2007-2015)")
> sdk + scale_color_manual(values = c("lightcoral", "purple"))
> #Bottom Left
> sdres <- plot(out.mbal15r2, type = "balance", stat = "sd",
>               main = "Covariate Balance\nMean\n(2008-2015)")
> sdres + scale_color_manual(values = c("lightcoral", "purple"))
> 
> #Appendix C.1.6--------------------------------------------------------------
> coefplot(M, drop = '(Intercept)',
>          pt.join = TRUE, ref = c("2009" = 3), ref.line = TRUE,
>          main = "Effect of Suffrage on Citizenship Acquisition\n2010 Enfranchised Origins in Spain",
>          xlab = "Year",
>          dict = c("year::2008"="2008", "year::2010"="2010",
>                   "year::2011"="2011", "year::2012"="2012"))
> modelsummary(M,
>              stars = T,
>              coef_map = c("2008","2009","2010","2011","2012"),
>              gof_omit = 'DF|Deviance|R2|Log.Lik.|BIC|Std. Errors')
> 
> 

========================================================================= 
Libraries 
========================================================================= 

========================================================================= 
Data Preparation 
========================================================================= 

rename: renamed 8 variables (Acquisitions, GDPpc, Unemployment, Democracy, Language, …)

NOTE: Log Print Time:  2024-03-16 13:59:57.981483 
NOTE: Elapsed Time: 2.56315875053406 mins 

========================================================================= 
Replicating Analysis 
========================================================================= 

Call:
tjbal.default(data = spain12, Y = "Acquisitions", D = "treat", 
    X = c("Colony", "Language", "Unemployment"), X.avg.time = list(c(2008), 
        c(2008), c(2008:2009)), index = c("scou", "year"), demean = T, 
    estimator = "mean", vce = "boot", conf.lvl = 0.9, nsims = 1000, 
    parallel = F, seed = 1013)

   ~ by Period (including Pre-treatment Periods):
      ATT   S.E. z-score CI.lower CI.upper p.value n.Treated
2008    0  130.8   0.000   -215.1    215.1  1.0000         7
2009    0  130.8   0.000   -215.1    215.1  1.0000         7
2010 4039 2234.5   1.808    363.4   7714.1  0.0707         7
2011 2141  884.3   2.421    686.1   3595.0  0.0155         7
2012 1637 1057.7   1.547   -103.0   3376.5  0.1218         7

Average Treatment Effect on the Treated:
      ATT S.E. z-score CI.lower CI.upper p.value
[1,] 2605 1010    2.58    944.5     4266  0.0099

NOTE: Log Print Time:  2024-03-16 14:03:40.804202 
NOTE: Elapsed Time: 3.71371198495229 mins 

Call:
tjbal.default(data = spain12, Y = "Acquisitions", D = "treat", 
    X = c("Colony", "Language", "Unemployment", "Democracy", 
        "GDPpc"), X.avg.time = list(c(2008), c(2008), c(2008:2009), 
        c(2008:2009), c(2008:2009)), index = c("scou", "year"), 
    demean = T, estimator = "mean", vce = "boot", conf.lvl = 0.9, 
    nsims = 1000, parallel = F, seed = 1013)

   ~ by Period (including Pre-treatment Periods):
      ATT   S.E. z-score CI.lower CI.upper p.value n.Treated
2008    0  134.3   0.000   -220.9    220.9  1.0000         7
2009    0  134.3   0.000   -220.9    220.9  1.0000         7
2010 4052 2217.0   1.828    405.3   7698.6  0.0676         7
2011 2159  849.0   2.543    762.6   3555.5  0.0110         7
2012 1679 1035.2   1.622    -23.5   3382.0  0.1048         7

Average Treatment Effect on the Treated:
      ATT  S.E. z-score CI.lower CI.upper p.value
[1,] 2630 977.3   2.691     1022     4238  0.0071

NOTE: Log Print Time:  2024-03-16 14:03:41.35377 
NOTE: Elapsed Time: 0.549567937850952 secs 

Call:
tjbal.default(data = spain7, Y = "Acquisitions", D = "treat", 
    X = c("Colony", "Language", "Unemployment", "Democracy", 
        "GDPpc"), X.avg.time = list(c(2008), c(2008), c(2007:2009), 
        c(2007:2009), c(2007:2009)), index = c("scou", "year"), 
    demean = T, estimator = "mean", vce = "boot", conf.lvl = 0.9, 
    nsims = 1000, parallel = F, seed = 1013)

   ~ by Period (including Pre-treatment Periods):
         ATT   S.E. z-score CI.lower CI.upper p.value n.Treated
2007  -731.7  366.1 -1.9986 -1333.91   -129.5  0.0456         7
2008   335.2  210.6  1.5919   -11.16    681.6  0.1114         7
2009   396.5  269.0  1.4741   -45.94    839.0  0.1405         7
2010  4472.9 2414.8  1.8523   500.86   8444.9  0.0640         7
2011  2600.2  957.9  2.7144  1024.55   4175.8  0.0066         7
2012  2120.5  969.1  2.1881   526.46   3714.5  0.0287         7
2013 10475.8 3358.3  3.1194  4951.95  15999.7  0.0018         7
2014 -1381.3 2432.2 -0.5679 -5381.99   2619.4  0.5701         7
2015 -2468.9 2723.2 -0.9066 -6948.19   2010.3  0.3646         7

Average Treatment Effect on the Treated:
      ATT S.E. z-score CI.lower CI.upper p.value
[1,] 2637 1052   2.506    906.3     4367  0.0122

NOTE: Log Print Time:  2024-03-16 14:03:41.87657 
NOTE: Elapsed Time: 0.522799968719482 secs 

Call:
tjbal.default(data = spain15, Y = "Acquisitions", D = "treat", 
    X = c("Perm.Res"), X.avg.time = list(c(2008:2009)), index = c("scou", 
        "year"), demean = T, estimator = "mean", vce = "boot", 
    conf.lvl = 0.9, nsims = 1000, parallel = F, seed = 1013)

   ~ by Period (including Pre-treatment Periods):
         ATT   S.E. z-score CI.lower CI.upper p.value n.Treated
2008  -120.4  131.2 -0.9177  -336.28    95.43  0.3588         7
2009   120.4  131.2  0.9177   -95.43   336.28  0.3588         7
2010  4026.6 2186.7  1.8414   429.79  7623.39  0.0656         7
2011  1752.3  834.9  2.0989   379.05  3125.57  0.0358         7
2012  1074.7 1235.0  0.8702  -956.70  3106.10  0.3842         7
2013  7163.1 4142.2  1.7293   349.78 13976.41  0.0838         7
2014 -2831.5 2968.5 -0.9539 -7714.25  2051.17  0.3401         7
2015 -3938.6 3276.0 -1.2023 -9327.12  1449.90  0.2293         7

Average Treatment Effect on the Treated:
      ATT S.E. z-score CI.lower CI.upper p.value
[1,] 1208 1539  0.7846    -1324     3740  0.4327

NOTE: Log Print Time:  2024-03-16 14:03:42.412463 
NOTE: Elapsed Time: 0.535892963409424 secs 

========================================================================= 
Log End Time: 2024-03-16 14:03:58.363146 
Log Elapsed Time: 0 00:06:34 
========================================================================= 
